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Spacecraft  Trajectory  Estimation 
Using  a  Sampled-Data  Extended  Kalman  Filter 
with  Range- Only  Measurements 


R.  Scott  Erwin  and  Dennis  S.  Bernstein 


Abstract — Determining  the  trajectory  of  a  body  orbiting 
the  Earth  is  a  fundamental  task  in  astrodynamics.  In  this 
paper  we  use  a  sampled-data  extended  Kalman  filter  to 
estimate  the  trajectory  of  a  target  satellite  when  only  range 
measurements  are  available  from  a  constellation  of  orbiting 
spacecraft.  We  consider  the  ability  of  the  filter  to  acquire 
the  target  satellite  under  time-sparse  measurements,  and  to 
estimate  the  eccentricity  and  inclination  of  the  target  satellite’s 
orbit.  Our  goal  is  to  quantify  tradeoffs  among  acquisition  time, 
tracking  accuracy,  and  measurement  sample  rate.  In  addition, 
when  the  orbits  of  the  observing  spacecraft  are  all  equatorial, 
it  is  found  that  inclination  maneuvers  of  the  target  satellite 
are  unobservable. 

1.  Introduction 

The  problem  of  estimating  the  full  state  of  a  dynamical 
system  based  on  limited  measurements  is  of  extreme  impor¬ 
tance  in  many  applications.  For  the  case  of  a  linear  system 
with  known  dynamics,  the  classical  Kalman  filter  provides 
an  optimal  solution  [1,2],  However,  state  estimation  for 
nonlinear  systems  remains  a  problem  of  intense  research 
interest. 

Besides  their  value  in  estimating  the  state  of  a  system 
with  nonlinear  dynamics,  nonlinear  estimators  can  also  be 
used  to  estimate  constant  states  that  represent  parameters. 
Consequently,  nonlinear  filters  are  useful  for  system  iden¬ 
tification  [3].  One  of  the  key  issues  that  arises  in  this 
application  is  parameter  bias,  a  longstanding  problem  [4]. 

Optimal  nonlinear  filters  have  been  studied  [5],  but  are 
often  infinite  dimensional  and  thus  are  difficult  to  imple¬ 
ment.  Within  a  deterministic  setting,  nonlinear  observers 
have  been  developed  for  systems  of  special  structure  [6, 
7],  Consequently,  except  for  systems  of  special  structure, 
approximate  filters  are  usually  implemented  in  practice. 

There  are  two  main  approaches  to  approximate  nonlinear 
filtering.  The  first  approach  is  based  on  a  linearization  of 
the  nonlinear  dynamics  and  measurement  mapping.  For 
example,  the  extended  Kalman  filter  uses  the  nonlinear 
dynamics  to  propagate  the  state  estimate  while  using  the 
linearized  dynamics  and  linearized  output  map  to  propagate 
the  pseudo-error  covariance.  The  extended  Kalman  filter  is 
often  highly  effective,  and  documented  applications  cover 
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an  extraordinarily  broad  range  of  disciplines,  from  motor 
control  to  weather  forecasting  [8,9], 

Although  the  extended  Kalman  filter  was  originally 
conceived  within  a  stochastic  setting,  recent  research  has 
provided  a  foundation  for  viewing  the  extended  Kalman 
filter  as  a  deterministic  observer  [10-12],  The  idea  is  to 
adopt  or  modify  the  formalisms  of  the  extended  Kalman 
filter  while  determining  conditions  that  ensure  stability  and 
convergence.  Although  the  sufficient  conditions  are  often 
conservative  for  specific  applications,  these  results  provide 
an  intellectual  rationalization  for  the  extended  Kalman  filter 
formalism 

The  second  approach  to  approximate  nonlinear  state 
estimation  foregoes  an  explicit  update  of  the  state  estimate 
error  covariance  in  favor  of  a  collection  of  filters  whose 
response  is  used  to  approximate  the  state  estimate  error 
covariance.  These  statistical  approaches  include  the  particle, 
unscented,  and  ensemble  Kalman  filters  [13-15]. 

The  present  paper  is  concerned  with  state  estimation  for 
satellite  trajectory  estimation,  which,  for  unforced  motion, 
is  equivalent  to  orbit  determination  [16],  Since  orbital 
dynamics  are  nonlinear,  nonlinear  estimation  techniques  are 
needed.  A  wide  variety  of  problems  can  be  considered 
based  on  the  type  of  data  that  are  available,  including  angle 
(azimuth  and  elevation),  range,  and  range  rate.  The  use  of 
angle-only  data  is  considered  in  [17],  which  develops  a 
specialized  filter  to  exploit  the  monotonicity  of  angles  in 
orbital  motion.  Issues  that  arise  in  the  use  of  range-rate 
(doppler)  measurements  are  discussed  in  [18, 19], 

Orbit  estimation  with  measurements  provided  by  a  con¬ 
stellation  of  satellites  is  considered  in  [20,21].  One  sce¬ 
nario  (TDRSS)  considers  the  use  of  observing  satellites  in 
circular,  equatorial,  geosynchronous  orbits  to  track  satellites 
in  low-Earth  orbit,  while  another  scenario  (GPS)  involves 
the  use  of  a  constellation  of  satellites  with  pseudo-range 
measurements  (range  measurements  with  clock  error  biases) 
to  determine  the  location  of  the  user. 

In  the  present  paper  we  consider  the  use  of  a  constellation 
of  satellites  in  low-Earth  orbit  to  track  a  satellite  in  geosyn¬ 
chronous  orbit.  Since  the  observing  satellites  have  much 
smaller  period  than  the  target  satellite  in  geosynchronous 
orbit,  we  must  account  for  blockage  by  the  Earth,  and 
thus  the  number  of  available  measurements  varies  with 
time.  We  are  particularly  interested  in  the  ability  of  the 
observing  satellites  to  acquire  and  track  the  target  satellite 


when  measurements  are  available  at  low  frequency,  that  is, 
with  a  large  sample  interval. 

Unlike  the  study  in  [20,21],  which  considers  either  angle 
(two  observations  per  satellite),  angle  and  range  (three 
observations  per  satellite),  and  angle,  range,  and  range  rate 
(four  observations  per  satellite),  we  consider  the  case  in 
which  only  range  measurements  are  available  (one  obser¬ 
vation  per  satellite).  In  addition,  the  sample  rate  in  [20]  was 
chosen  to  be  33  Hz,  while  we  are  interested  in  the  ability 
to  track  under  time-sparse  measurements,  available  on  a 
scale  of  only  minutes  or  perhaps  hours.  This  constraint  is 
motivated  by  the  need  for  satellites  to  simultaneously  track 
a  large  number  of  objects. 

As  in  [20,21],  we  employ  the  sampled-data  (continuous- 
discrete)  extended  Kalman  filter  [2,  p.  188],  This  extended 
Kalman  filter  involves  continuous-time  propagation  of  the 
state  estimate  as  well  as  the  pseudo-error  covariance  be¬ 
tween  measurements  and  data  updates.  In  practice,  the 
continuous-time  state  and  covariance  propagation  can  be 
implemented  online  with  high-resolution  integration  to  ac¬ 
curately  follow  the  nonlinear  dynamics. 

2.  Equations  of  Motion 

We  consider  a  single  body,  called  the  target,  orbiting 
the  Earth.  Throughout  this  study  we  assume  that  the  Earth 
is  spherical  and  uniform  with  an  ideal  gravitational  field. 
Except  for  possible  thrusting  by  the  target  itself,  we  ignore 
all  perturbing  forces  such  as  drag.  The  position  vector  r  of 
the  target  with  respect  to  the  center  of  the  Earth  satisfies 

r  =  — ^  r  +  w,  (1) 


where  r  =  |  r  |  is  the  distance  from  the  satellite  to  the  center 
of  the  Earth,  w  denotes  forces  due  to  thrusting  per  unit  mass 
acting  on  the  target,  and  /i  =  398, 600  km3/s2  is  the  Earth’s 
gravitational  parameter.  The  specific  thrust  w  is  zero  unless 

the  target  is  actively  maneuvering.  Introducing  the  velocity 
-»  A  -4  ... 

vector  v  =  r,  we  can  rewrite  (1)  as 

~r  =  v,  (2) 

v  =  —3-  r  +  w.  (3) 

To  cast  the  dynamics  (3)  in  terms  of  coordinates,  we 
introduce  an  inertial  reference  frame  I.  It  is  traditional  to 
choose  the  inertial  reference  frame  so  that  the  X-axis  points 
toward  the  Sun  on  the  first  day  of  spring  (the  vernal  equinox 
line),  the  Z-axis  points  through  the  North  pole  of  the  Earth 
along  the  axis  of  rotation,  and  the  Y -axis  completes  a  right- 
handed  coordinate  system.  This  description  is  approximate 
since  the  Earth’s  rotational  axis  is  not  fixed  inertially  and 
since  the  stars  move  inertially  as  well  [22,  pp.  150-153], 
However,  such  details  do  not  play  a  role  in  the  subsequent 
analysis. 


Resolving  r,v,  and  w  in  I  according  to 
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Note  that  r  =  \J  x1  +  y2  +  z 2.  We  can  rewrite  (4)  as 
X  =  /(X)  +  w, 


where 
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The  vector  X  provides  a  complete  representation  of 
the  target’s  state,  which  is  characterized  by  the  position 
and  velocity.  When  the  satellite  is  moving  along  an  orbit, 
such  as  a  circle,  ellipse,  parabola,  or  hyperbola,  it  is  often 
useful  to  represent  the  satellite  motion  in  terms  of  the  6 
orbital  parameters  given  by  the  perigee  distance  rp,  the 
eccentricity  e,  the  right  ascension  of  the  ascending  node 
f2,  the  inclination  i,  the  argument  of  periapsis  at,  and  the 
true  anomaly  v.  The  orbital  elements  rp  and  e  fix  the  shape 
of  the  orbit,  while  the  angles  f 2,  i,  and  oj  comprise  a  (3, 1,  3) 
sequence  of  Euler  rotations  that  transform  the  inertial  frame 
to  the  inertially  fixed  frame.  The  true  anomaly  v(t)  is  a 
time-dependent  parameter  that  keeps  track  of  the  position 
of  the  satellite  along  its  orbit.  The  nonlinear  transformations 
that  convert  position  and  velocity  into  orbital  elements  and 
vice  versa  are  given  in  [22]. 


3.  Measurement  Model 


For  trajectory  estimation,  we  assume  that  range  measure¬ 
ments  are  available  from  p  satellites  at  times  t  =  kh,  where 
k  =  1,2, ... .  Letting  Xi,  yi,  Zi  denote  the  inertial  coordi¬ 
nates  of  the  ith  satellite,  assumed  to  be  known  accurately, 
the  measurement  Y  =  Y(kh)  £  Rp  is  given  by  (omitting 
the  argument  kh  on  the  right-hand  side) 


Y(kh) 


di(x,y,z,x1,y1,z1) 


dp(x,y,z,xp,yp,Zp) 


+  v, 


(7) 


where,  for  i  =  1, . . .  ,p, 


A 


di(x,y,z,Xi,yi,Zi)  = 

[{x  -  Xi)2  +  (y-  yi)2  +  (z-  Zi)2}1/2  (8) 


is  the  distance  from  the  ith  satellite  to  the  target,  and  v  €  Rp 
denotes  measurement  noise. 

The  measurement  di  is  assumed  to  be  unavailable  when 
the  line-of-sight  path  between  the  ?'th  satellite  and  the  target 
is  blocked  by  the  Earth.  To  determine  blockage,  we  note  that 
the  Earth’s  surface  blocks  the  path  from  the  ith  satellite  to 
the  target  if  and  only  if  there  exists  a  £  [0, 1]  such  that 

Di(a)  <  Re ,  where  Re  =  6378  km  is  the  radius  of  the 
Earth  and 

Di{a)  =  \/  (xi  +  ax)2  +  (t/*  +  ay)2  +  +  az)2. 

The  smallest  value  of  Di{a)  as  attained  for  a  =  an,  where 

A  XiX  +  UiU  +  ZiZ 

O^i  o  i  o  i  o 

x-  +  jr  +  zz 

Hence,  we  compute  on,  ascertain  whether  a*  lies  in  the 
interval  [0,1],  and  then  check  the  blockage  condition 

Di(ai)  <  Re- 

4.  Sampled-Data  Extended  Kalman  Filter 

Since  the  equations  of  motion  (1)  are  nonlinear,  we 
consider  an  extended  Kalman  filter.  In  addition,  since  we 
assume  that  measurements  Y  are  available  with  a  specified 
sample  interval  of  h  sec,  we  consider  a  sampled-data 
extended  Kalman  filter  with  data  update  performed  at  each 
time  t  =  kh.  The  sampled-data  (“continuous-discrete”) 
extended  Kalman  filter  is  given  in  [2,  p.  188].  We  use 
”  and  “+”  to  denote  state  estimates  before  and  after  data 
updates,  respectively. 

4.1.  Forecast  Step 

The  forecast  (data-free)  step  of  the  sampled-data  ex¬ 
tended  Kalman  filter  consists  of  the  state-estimate  propa- 


gation 

i(£)  =  f(X(t)),  £  G  [(fe  -  l)h,kh\, 

(9) 

where 

X  =  [  x  y  z  vx  vy  i>z  ] T , 

(10) 

as  well  as  the  pseudo-error  covariance  propagation 

P(t)  =  A(t)P(t)  +  P(£)iT(£)  +  Q,  £  e  [(fe  —  1  )h,  kh], 

(11) 


Equations  (9)  and  (11)  are  numerically  integrated  in  real 
time.  Since  there  is  no  data  injection  during  the  time  interval 
[(fe  —  1  )h,kh],  variable-step  integration  can  be  used  for 
efficiency  and  accuracy.  Let  X(kh—)  and  P(kh-)  denote 
the  values  of  X  and  P  at  the  right-hand  endpoint  of  the 
interval  [(fe  —  1  )h,kh\.  At  the  start  of  the  next  interval 
[kh,  ( k  +  1  )h],  the  initial  values  X(kh+)  and  P(kh+)  are 
determined  by  the  data  update  step.  The  overall  system  can 
be  viewed  as  a  sampled-data  system  in  which  continuous¬ 
time  dynamics  are  interrupted  by  instantaneous  state  jumps 
[23], 

4.2.  Data  Update  Step 

For  the  data  update  step,  the  linearized  measurement  map 
is  given  by 


Ck  = 


x(kh—)—xi(kh)  y(kh—)  —  yi(kh)  z{kh—)—z\(kh)  z-v 

di(fc)  di(k)  di(k) 

x(kh—)—xp(kh)  y(kh—)  —  yp(kh)  z(kh—)  —  zp(kh)  q 


1x3 


dp(k)  dp(k) 

where,  for  i  =  1, . . . 


dp  (fe) 


1x3 


di{k) 


A 


di(x(kh—),y(kh—),z(kh-),Xi(kh),yi(kh),  Zi(kh)). 


Furthermore,  the  data  update  gain  Kk  is  given  by 
Kk  =  P(kh-)C^[CkP(kh-)C^  +  R]~\ 
while  the  state-estimate  data  update  is  given  by 

X(kh+)  =  X{kh—)  +  Kk[Y(kh)  -  d(k)], 

where 


d(fe)  = 


dt(fe) 

dp(k)  J 


(12) 

(13) 


Finally, 

P(kh+)  =  (I  -  KkCk)P(kh—). 

The  values  X(kh+)  and  P(kh+)  are  used  to  initialize  (9), 
(11)  in  the  next  interval  [kh,  ( k  +  l)fe]. 


5.  Numerical  Examples 


where  A(t)  =  f'(X(t))  is  the  Jacobian  of  /  evaluated  along 
the  trajectory  of  (9).  The  Jacobian  A(t)  is  given  by 


Ait) 


03x3  h 
A)(£)  03x3 


where  (omitting  the  argument  t ) 


A0{t)  =  n 


3& _ ^  3  xy  3  xz 

3xy  3  y^ _ 3  yz 

3 xz  3yz  3z2  1 


where  r  =  \J x2  +  y2  +  z2. 


We  consider  the  case  in  which  satellites  in  low-Earth 
orbit  (LEO)  at  a  radius  of  6600  km  are  observing  a  target 
satellite  in  an  equatorial  geosynchronous  orbit  at  a  radius  of 
42,164  km.  We  assume  that  the  LEO  satellites  are  spaced 
uniformly  around  the  Earth  in  an  equatorial  orbit.  With 
this  arrangement,  4  is  the  smallest  number  of  satellites  for 
which  at  least  2  satellites  are  always  able  to  simultaneously 
view  the  target.  Target  tracking  with  as  few  as  2  satellites 
separated  by  a  true  anomaly  of  less  than  180  degrees  is 
also  possible  as  long  as  measurements  are  guaranteed  to 
be  available  when  the  target  is  in  the  field  of  view  of 
both  satellites.  However,  for  simplicity,  we  assume  the 
availability  of  4  uniformly  spaced  LEO  satellites  with  range 


measurements  available  (subject  to  blocking  by  the  Earth)  at 
a  fixed  sample  interval  of  ft  sec.  All  satellite  measurements 
(blocked  or  not)  are  assumed  to  occur  simultaneously. 

Assuming  perfect  knowledge  of  the  initial  condition  and 
assuming  that  the  target  is  not  maneuvering,  it  is  possible 
to  track  the  target  with  arbitrary  accuracy  without  the  use 
of  measurements.  In  practice,  perturbing  forces  such  as 
drag  must  also  be  estimated;  however,  these  forces  are 
not  considered  in  this  study.  When  the  initial  state  is 
unknown  or  when  the  target  is  maneuvering,  measurements 
are  needed  to  track  the  target.  We  consider  these  cases 
separately. 

6.  Target  Acquisition 

We  first  consider  the  ability  of  the  sampled-data  Kalman 
filter  to  acquire  the  target,  that  is,  to  locate  the  target  despite 
initial  position  errors.  In  all  cases  we  set  Q  =  [  o  /3  ]  ,  and 

=  o. 

First,  we  set  the  sample  interval  to  be  ft,  =  1  sec,  and  we 
consider  initial  estimates  that  are  erroneous  by  1  degree  and 
110  degrees.  We  assume  perfect  (nonnoisy)  measurements 
and  set  R  =  0.01/  in  (12).  Figure  1  compares  the  perfor¬ 
mance  of  the  filter  for  both  initial  estimates.  The  ultimate 
tracking  accuracy  in  both  cases  is  determined  by  numerical 
resolution  in  computing  the  state  estimates.  Convergence  of 
the  filter  is  not  global;  in  fact,  the  filter  fails  to  converge  for 
initial  true  anomaly  errors  greater  than  about  120  degrees. 

Next,  we  introduce  gaussian  measurement  noise  with  a 
standard  deviation  of  0.1  km,  which  corresponds  to  R  = 
0.01/  in  (12).  For  initial  estimates  that  are  erroneous  by  1 
degree  and  110  degrees.  Figure  2  shows  that  the  position- 
estimate  error  reaches  a  level  that  is  consistent  with  the 
measurement  accuracy. 

For  an  initial  true  anomaly  estimate  that  is  erroneous  by 
110  degrees,  the  position  estimates  are  shown  in  Figure  3. 
The  estimator  approaches  the  vicinity  of  the  target  within 
about  10  sec. 

Next  we  consider  the  ability  of  the  filter  to  acquire  the 
target  under  time-sparse  measurements  with  a  measurement 
standard  deviation  of  0.1  km.  For  an  initial  true  anomaly 
error  of  10  degrees.  Figure  4  shows  the  position-estimate 
errors  for  ft  =  1, 10,  50, 100  sec.  In  each  case,  the  estimator 
acquires  the  target  in  about  10  data  assimilation  steps,  with 
ultimate  accuracy  independent  of  the  sample  interval. 

7.  Intermeasurement  Tracking  Accuracy 

Next,  we  assess  the  ability  of  the  filter  to  track  the  target 
along  its  orbit.  To  see  how  the  position  estimate  degrades 
between  data  updates,  we  consider  an  initial  true  anomaly 
error  of  10  degrees  and  a  sample  interval  of  50  sec  with 
measurement  noise  having  a  standard  deviation  of  0.1  km. 
Figure  5  shows  the  growth  of  the  position  error  between 
measurements  as  well  as  the  position-error  reduction  that 
occurs  due  to  data  injection. 
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Fig.  1.  Target  position-  and  velocity-estimate  errors  with  initial 
true  anomaly  errors  of  1  degree  and  110  degrees.  The  range  data 
are  measured  with  a  sample  interval  h  =  1  sec  from  4  LEO 
satellites,  and  perfect  (nonnoisy)  measurements  are  assumed  with 
R  =  0.01 1  in  the  filter  gain  expression  (12).  In  both  cases  the 
estimator  accurately  locates  the  target  within  about  50  sec. 
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Fig.  2.  Target  position-  and  velocity-estimate  errors  with  initial 
true  anomaly  errors  of  1  degree  and  110  degrees.  The  range  data 
are  measured  with  a  sample  interval  h  =  1  sec  from  4  LEO 
satellites,  and  with  gaussian  measurement  noise  whose  standard 
deviation  is  0.1  km  and  thus  with  R  =  0.01/  in  the  filter  gain 
expression  (12).  In  both  cases  the  estimator  accuracy  corresponds 
to  the  measurement  error  level. 


8.  Eccentricity  Estimation 

We  now  consider  the  case  in  which  the  target  performs 
an  unknown  thrust  maneuver  that  changes  the  eccentricity 
of  its  orbit.  The  initial  true  anomaly  error  in  all  cases  is  10 
degrees.  In  particular,  the  target  is  initially  in  a  circular  orbit 
as  in  the  previous  examples.  At  time  t  =  100  sec,  the  target 
performs  a  1 -second  burn  that  produces  a  specific  thrust 
w  =  [0  .5  0]T  km/s2,  and,  at  time  t  =  200  sec,  the  target 
performs  a  1 -second  burn  that  produces  a  specific  thrust 
w  =  [0  .3  0]T  km/s2.  With  an  initial  eccentricity  of  e  =  0, 
corresponding  to  the  initial  circular  orbit,  the  eccentricity 
after  the  first  burn  is  e  «  .35,  while  the  eccentricity  after 
the  second  burn  is  e  ~  .59.  Assuming  measurement  noise 
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Fig.  4.  Target  position-estimate  error  for  sample  intervals  h  = 
1,20,50,100  with  range  data  measured  from  4  LEO  satellites 
with  a  standard  deviation  of  0.1  km.  In  each  case,  the  estimator 
acquires  the  target  in  about  10  data  assimilation  steps,  with  ultimate 
accuracy  independent  of  the  sample  interval. 
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Fig.  3.  Target  position  estimates  ’ x’  with  an  initial  true  anomaly 
error  of  110  degrees.  The  initial  location  of  the  target  is  at  3:00. 
The  range  data  are  measured  with  a  sample  interval  h  =  1  sec 
from  4  LEO  satellites  (whose  tracks  are  shown),  and  with  gaussian 
measurement  noise  whose  standard  deviation  is  0.1  km  and  thus 
with  R  =  0.01/  in  the  filter  gain  expression  (12).  The  estimator 
approaches  the  vicinity  of  the  target  within  about  10  sec.  The  Earth 
and  all  orbit  locations  are  drawn  to  scale. 


with  a  standard  deviation  of  0.01  km  and  measurements 
available  with  a  sample  interval  of  ft  =  1  sec,  the  estimated 
eccentricity  based  on  the  data  update  estimates  is  shown  in 
Figure  6.  The  same  scenario  is  repeated  with  ft  increased 
to  10  sec,  with  the  results  shown  in  Figure  7. 

9.  Inclination  Estimation 

We  now  consider  the  case  in  which  the  target  performs 
an  unknown  thrust  maneuver  that  changes  its  inclination. 
The  initial  true  anomaly  error  in  all  cases  except  where 
noted  is  30  degrees.  In  particular,  the  target  is  initially  in 
a  circular  equatorial  orbit  as  in  the  previous  examples.  At 
time  t  =  100  sec,  the  target  performs  a  1 -second  burn  that 
produces  a  specific  thrust  w  =  [0  0  .5]T  km/s2,  and,  at 
time  t  =  200  sec,  the  target  performs  a  1 -second  burn  that 
produces  a  specific  thrust  w  =  [0  0  —  .2]T  km/s2.  With  an 
initial  inclination  of  i  =  0  rad,  corresponding  to  the  initial 
equatorial  orbit,  the  inclination  after  the  first  burn  is  t  « 
0.16  rad,  while  the  inclination  after  the  second  burn  is  i  ~ 
0.097  rad.  Assuming  measurement  noise  with  a  standard 
deviation  of  0.01  km,  the  estimated  inclination  based  on 
the  data  update  estimates  is  shown  in  Figure  8. 

Figure  8  shows  that,  after  a  transient,  the  filter  correctly 
converges  to  the  target’s  inclination  of  0  rad.  However,  the 
filter  fails  to  detect  the  changes  in  inclination  due  to  the 
target’s  maneuvers.  This  failure  suggests  that  the  out-plane- 
maneuver  is  due  to  a  lack  of  observability  by  the  observing 
satellites.  In  fact,  the  numerical  rank  of  the  observability 
matrix  formed  from  (A(kh—),Ck)  is  found  to  be  4. 
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Fig.  5.  Target  position-estimate  error  with  an  initial  true  anomaly 
error  of  10  degrees,  a  sample  interval  of  h  =  50  sec,  and  a 
measurement  noise  standard  deviation  of  0.1  km.  The  growth  of 
the  position  error  between  measurements  can  be  seen,  as  well  as 
the  position-error  reduction  that  occurs  due  to  data  injection. 


Next,  we  slightly  change  the  orbit  of  the  first  observing 
satellite  by  giving  it  an  inclination  of  —0.1  rad.  Figure 
9  shows  a  strong  transient  that  prevents  the  filter  from 
estimating  the  initial  0  rad  inclination  despite  the  fact  that 
the  initial  inclination  estimate  is  correct,  followed  by  a  noisy 
estimate  of  the  inclination  after  the  first  burn,  followed, 
finally,  by  a  biased  estimate  of  the  inclination  after  the 
second  burn. 

Next,  we  also  change  the  orbit  of  the  second  observing 
satellite  by  giving  it  an  inclination  —0.2  rad.  After  an  initial 
transient.  Figure  10  shows  an  improved  ability  to  estimate 
the  true  inclination. 

Finally,  we  increase  the  sample  interval  to  ft  =  10  sec. 
In  this  case,  the  filter  diverges.  However,  for  an  initial  true 
anomaly  error  of  5  degrees,  it  can  be  seen  from  Figure  1 1 
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Fig.  6.  Estimated  eccentricity.  The  target  performs  unknown 
1 -second  bums  at  t  =  100  sec  and  t  =  200  sec.  The  initial 
eccentricity  is  e  =  0,  corresponding  to  the  initial  circular  orbit, 
while  the  eccentricity  after  the  first  burn  is  e  ~  .35,  and  the 
eccentricity  after  the  second  burn  is  e  ~  .59.  Assuming  a  1 -second 
measurement  interval  and  measurement  noise  with  a  standard 
deviation  of  0.01  km,  the  estimated  eccentricity  follows  the  true 
eccentricity.  The  apparent  bias  for  the  initial  eccentricity  e  =  0  is 
an  artifact  of  the  constraint  e  >  0.  The  full  extent  of  the  initial 
estimate  transient  is  not  shown. 
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Fig.  7.  Estimated  eccentricity.  This  simulation  is  analogous  to 
Figure  6,  where  now  the  sample  interval  is  h  =  10  sec.  Again,  the 
full  extent  of  the  initial  estimate  transient  is  not  shown. 


Fig.  8.  Estimated  inclination.  Although  the  filter  correctly  con¬ 
verges  to  the  target’s  initial  inclination  of  0  rad,  the  filter  fails  to 
detect  the  changes  in  inclination  due  to  the  target’s  maneuvers. 
Numerical  tests  suggest  a  lack  of  observability  by  the  observing 
satellites. 
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Fig.  9.  Estimated  inclination.  In  this  case,  the  orbit  of  the 
first  observing  satellite  by  given  an  inclination  of  —0.1  rad.  A 
strong  transient  prevents  the  filter  from  estimating  the  initial  0  rad 
inclination  despite  the  fact  that  the  initial  inclination  estimate  is 
correct.  The  estimate  of  the  inclination  after  the  first  maneuver  is 
noisy,  while  the  estimate  of  the  inclination  after  the  second  burn 
is  biased. 


that  the  filter  can  detect  changes  in  the  target’s  inclination. 

10.  Concluding  Remarks 

Under  idealized  assumptions  on  the  astrodynamics  of 
bodies  orbiting  the  Earth,  we  developed  a  sampled-data 
Kalman  filter  for  range-only  observations  provided  by  a 
constellation  of  4  low-Earth  orbiting  satellites  in  circular, 
equatorial,  geosynchronous  orbits.  We  studied  the  ability 
of  the  filter  to  acquire  and  track  a  target  satellite  in 
geosynchronous  orbit  as  a  function  of  the  sample  interval, 
initial  true  anomaly  error,  and  measurement  noise  standard 
deviation. 

This  study  complements  previous  studies  that  have  con¬ 
sidered  combinations  of  angle,  range,  and  range-rate  mea¬ 


surements  by  considering  range  measurements  alone  and  by 
considering  the  effects  of  infrequent  measurements. 

A  surprising  discovery  of  our  study  is  the  apparent 
inability  of  a  constellation  of  4  LEO  satellites  to  track 
a  maneuvering  target  satellite  that  changes  its  inclination. 
Naive  numerical  tests  based  on  linear  time-invariant  no¬ 
tions  suggest  that  the  dynamics  of  the  target  satellite  are 
unobservable  using  the  available  measurements.  Therefore, 
future  research  will  seek  to  apply  nonlinear  observability 
tests  [24-26]  to  better  understand  how  observability  of  the 
target’s  dynamics  depends  on  the  geometry  of  the  observing 
constellation  and  the  types  of  available  measurements. 
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Fig.  10.  Estimated  inclination.  In  addition  to  the  inclination  of  the 
orbit  of  the  first  observing  satellite,  the  second  observing  satellite 
is  given  an  inclination  of  —0.2  rad.  After  an  initial  transient,  the 
filter  provides  improved  estimates  of  the  target’s  inclination. 
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Fig.  11.  Estimated  inclination.  In  addition  to  the  inclinations  of 
the  orbits  of  the  first  and  second  observing  satellites,  the  sample 
interval  is  increased  to  h  =  10  sec.  Despite  this  constraint,  the 
filter  retains  the  ability  to  detect  changes  in  the  target’s  inclination. 
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